list.of.packages <- c("R.matlab")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(R.matlab)

# FUNCTION: EIGENVALUE-CORRECTION ######################################
EV_Correction <- function( A ){
 EV = eigen(A)
 Eval = EV$values
 Eval = (Eval >= 0) * Eval + (Eval < 0 ) * 0
 Evec = EV$vectors
 return( (Evec) %*% diag(Eval) %*% t(Evec) )
}

first_period = 427
last_period = 546
period = first_period:last_period
calendar = function(u){paste(floor((u+5)/12)+1964,"/",(u+5)%%12+1,sep="")}
paste(calendar(first_period)," - ",calendar(last_period),sep="")
##################################################################
# FAMA-FRENCH 4 FACTORS (BEGIN)
##################################################################
Data_factors <- readMat("Wspace_Fact.mat")$F
MKT = Data_factors[,1] # 3-Factor
SMB = Data_factors[,2] # 3-Factor
HML = Data_factors[,3] # 3-Factor
MOM = Data_factors[,4] # 4-Factor
##################################################################
# FAMA-FRENCH 4 FACTORS (ENDED)
##################################################################

##################################################################
# RETURNS (BEGIN)
##################################################################
R_25FF <- readMat("Wspace_25FF.mat")$R        # 25 Fama-French Portfolios
R_44IN <- readMat("Wspace_44Indu.mat")$R      # 44 Industry Portfolios
R_ALL <- readMat("Wspace_CRSPCMST_ret.mat")$R # 9936 Stocks Matching CRSP & Compustat
R_RF <- readMat("RiskFree.mat")$Rf            # Risk-Free Rate
##################################################################
# RETURNS (ENDED)
##################################################################

##################################################################
# GET ONLY THE BALANCED PORTION OF ALL THE STOCKS (BEGIN)
##################################################################
R_25FF = R_25FF[period,]
R_44IN = R_44IN[period,]
list_balanced = NULL
for(idx in 1:(dim(R_ALL)[2])){
 if(sum(R_ALL[period,idx]==0)<1){
  list_balanced = c(list_balanced,idx) 
 }
}
R_ALL = R_ALL[period,list_balanced]
##################################################################
# GET ONLY THE BALANCED PORTION OF ALL THE STOCKS (ENDED)
##################################################################

##################################################################
# CONSTRUCT X & Y (BEGIN)
##################################################################
R = R_ALL # THIS LINE CAN BE CHANGED TO R_25FF, R_44IN OR R_ALL ###################################################################
T = dim(R)[1]
N = dim(R)[2]
Y = R - matrix(R_RF[period],T,N)
X1 = matrix(MKT[period],T,N)
X2 = matrix(SMB[period],T,N)
X3 = matrix(HML[period],T,N)
##################################################################
# CONSTRUCT X & Y (ENDED)
##################################################################

##################################################################
# WITHIN-TRANSFORMATION (BEGIN)
##################################################################
Y = t((Y - t(matrix(colMeans(Y),N,T)))[-1,])
X1 = t((X1 - t(matrix(colMeans(X1),N,T)))[-1,])
X2 = t((X2 - t(matrix(colMeans(X2),N,T)))[-1,])
X3 = t((X3 - t(matrix(colMeans(X3),N,T)))[-1,])
T = T - 1
##################################################################
# WITHIN-TRANSFORMATION (ENDED)
##################################################################

##################################################################
# OLS (BEGIN)
##################################################################
vecY = c(Y)
vecX = cbind(c(X1),c(X2),c(X3))
beta_hat = solve(t(vecX)%*%vecX)%*%t(vecX)%*%vecY
beta_hat
##################################################################
# OLS (ENDED)
##################################################################

##################################################################
# RESIDUALS (BEGIN)
##################################################################
U_hat = Y - beta_hat[1,1]*X1 - beta_hat[2,1]*X2 - beta_hat[3,1]*X3
##################################################################
# RESIDUALS (ENDED)
##################################################################



########################################################################
# ESTIMATE IID OMEGA (BEGIN)
########################################################################
Omega_hat_iid = matrix(0,3,3)
Omega_hat_iid[1,1] = t(matrix(c(X1*U_hat),,1)) %*% matrix(c(X1*U_hat),,1)
Omega_hat_iid[2,1] = t(matrix(c(X2*U_hat),,1)) %*% matrix(c(X1*U_hat),,1)
Omega_hat_iid[3,1] = t(matrix(c(X3*U_hat),,1)) %*% matrix(c(X1*U_hat),,1)
Omega_hat_iid[1,2] = t(matrix(c(X1*U_hat),,1)) %*% matrix(c(X2*U_hat),,1)
Omega_hat_iid[2,2] = t(matrix(c(X2*U_hat),,1)) %*% matrix(c(X2*U_hat),,1)
Omega_hat_iid[3,2] = t(matrix(c(X3*U_hat),,1)) %*% matrix(c(X2*U_hat),,1)
Omega_hat_iid[1,3] = t(matrix(c(X1*U_hat),,1)) %*% matrix(c(X3*U_hat),,1)
Omega_hat_iid[2,3] = t(matrix(c(X2*U_hat),,1)) %*% matrix(c(X3*U_hat),,1)
Omega_hat_iid[3,3] = t(matrix(c(X3*U_hat),,1)) %*% matrix(c(X3*U_hat),,1)
Omega_hat_iid = 1/(N*T) * Omega_hat_iid
########################################################################
# ESTIMATE IID OMEGA (ENDED)
########################################################################



########################################################################
# MENZEL'S BOOTSTRAP (BEGIN)
########################################################################

set.seed(0)

# MENZEL TUNING PARAMETER ##############################################
num_boot = 2000

# MAMEN RANDOM VARIABLE ################################################
mammen <- function(para){ a = para[1]; b = para[2]; p = para[3]; return( (p*a+(1-p)*b)^2 + (p*a^2+(1-p)*b^2-1)^2 + (p*a^3+(1-p)*b^3-1)^2 ); }
mammen_para = nlm(mammen,c(-1,1,0.5))$estimate
rmammen <- function(n,para){(para[1]-para[2])*(runif(n)<para[3])+para[2]}

X1Ubar = mean(X1*U_hat)
X2Ubar = mean(X2*U_hat)
X3Ubar = mean(X3*U_hat)
a1 = rowMeans(X1*U_hat)
a2 = rowMeans(X2*U_hat)
a3 = rowMeans(X3*U_hat)
g1 = colMeans(X1*U_hat)
g2 = colMeans(X2*U_hat)
g3 = colMeans(X3*U_hat)
w1 = X1*U_hat - matrix(a1,N,T) - t(matrix(g1,T,N))
w2 = X2*U_hat - matrix(a2,N,T) - t(matrix(g2,T,N))
w3 = X3*U_hat - matrix(a3,N,T) - t(matrix(g3,T,N))

S2_a1 = rowSums((X1*U_hat - X1Ubar)^2)/(N-1)
S2_a2 = rowSums((X2*U_hat - X2Ubar)^2)/(N-1)
S2_a3 = rowSums((X3*U_hat - X3Ubar)^2)/(N-1)
S2_g1 = colSums((X1*U_hat - X1Ubar)^2)/(T-1)
S2_g2 = colSums((X2*U_hat - X2Ubar)^2)/(T-1)
S2_g3 = colSums((X3*U_hat - X3Ubar)^2)/(T-1)
S2_w1 = sum((X1*U_hat - matrix(a1,N,T) - t(matrix(g1,T,N)) - X1Ubar)^2) / (N*T - N - T)
S2_w2 = sum((X2*U_hat - matrix(a2,N,T) - t(matrix(g2,T,N)) - X2Ubar)^2) / (N*T - N - T)
S2_w3 = sum((X3*U_hat - matrix(a3,N,T) - t(matrix(g3,T,N)) - X3Ubar)^2) / (N*T - N - T)

sigma2_a1 = max(c(0,S2_a1 - S2_w1/T))
sigma2_a2 = max(c(0,S2_a2 - S2_w2/T))
sigma2_a3 = max(c(0,S2_a3 - S2_w3/T))
sigma2_g1 = max(c(0,S2_g1 - S2_w1/N))
sigma2_g2 = max(c(0,S2_g2 - S2_w2/N))
sigma2_g3 = max(c(0,S2_g3 - S2_w3/N))
sigma2_w1 = S2_w1
sigma2_w2 = S2_w2
sigma2_w3 = S2_w3

eps_hat = U_hat - matrix(a1,N,T) - t(matrix(g1,T,N)) - mean(U_hat)
mu4_e = sqrt(2)*max(c(0.1,sqrt(mean(eps_hat^4)-mean(eps_hat)^4)))
kappa_a = mu4_e * log(T)
kappa_g = mu4_e * log(N)

D_a1 = 1*(T*sigma2_a1 >= kappa_a)
D_a2 = 1*(T*sigma2_a2 >= kappa_a)
D_a3 = 1*(T*sigma2_a3 >= kappa_a)
D_g1 = 1*(N*sigma2_g1 >= kappa_g)
D_g2 = 1*(N*sigma2_g2 >= kappa_g)
D_g3 = 1*(N*sigma2_g3 >= kappa_g)

lambda_a1 = (D_a1 * T * sigma2_a1) / (D_a1*T*sigma2_a1 + sigma2_w1)
lambda_a2 = (D_a2 * T * sigma2_a2) / (D_a2*T*sigma2_a2 + sigma2_w2)
lambda_a3 = (D_a3 * T * sigma2_a3) / (D_a3*T*sigma2_a3 + sigma2_w3)
lambda_g1 = (D_g1 * N * sigma2_g1) / (D_g1*N*sigma2_g1 + sigma2_w1)
lambda_g2 = (D_g2 * N * sigma2_g2) / (D_g2*N*sigma2_g2 + sigma2_w2)
lambda_g3 = (D_g3 * N * sigma2_g3) / (D_g3*N*sigma2_g3 + sigma2_w3)

boot_beta = NULL
for ( boot_iter in 1:num_boot ){

# DRAW BOOTSTRAP SAMPLES ###############################################
ki = sample(1:N,N,replace=TRUE)
kt = sample(1:T,T,replace=TRUE)

# STEP (a) #############################################################
a1_star = a1[ki]
a2_star = a2[ki]
a3_star = a3[ki]
g1_star = g1[kt]
g2_star = g2[kt]
g3_star = g3[kt]
# STEP (b) #############################################################
omega1_star = matrix(rmammen(N,mammen_para),N,T)
omega2_star = t(matrix(rmammen(T,mammen_para),T,N))
w1_star = omega1_star*omega2_star*w1[ki,kt]
w2_star = omega1_star*omega2_star*w2[ki,kt]
w3_star = omega1_star*omega2_star*w3[ki,kt]
# STEP (c) #############################################################
z1_star = lambda_a1^0.5*matrix(a1_star,N,T) + lambda_g1^0.5*t(matrix(g1_star,T,N)) + w1_star
z2_star = lambda_a2^0.5*matrix(a2_star,N,T) + lambda_g2^0.5*t(matrix(g2_star,T,N)) + w2_star
z3_star = lambda_a3^0.5*matrix(a3_star,N,T) + lambda_g3^0.5*t(matrix(g3_star,T,N)) + w3_star
# STEP (d) #############################################################
boot_beta = cbind(boot_beta, beta_hat + N*T*solve( t(cbind(matrix(c(X1),,1),matrix(c(X2),,1),matrix(c(X3),,1))) %*% cbind(matrix(c(X1),,1),matrix(c(X2),,1),matrix(c(X3),,1)) ) %*% matrix(c(mean(z1_star),mean(z2_star),mean(z3_star)),3,1))
}

SE_MEN = apply(boot_beta,1,sd)

########################################################################
# MENZEL'S BOOTSTRAP (ENDED)
########################################################################



########################################################################
# ESTIMATE PANEL NEWEY-WEST OMEGA (BEGIN)
########################################################################
Omega_hat_PNW = matrix(0,3,3)
Omega_hat_PNW[1,1] = t(matrix(c(X1*U_hat))) %*% matrix(c(X1*U_hat))
Omega_hat_PNW[2,1] = t(matrix(c(X2*U_hat))) %*% matrix(c(X1*U_hat))
Omega_hat_PNW[3,1] = t(matrix(c(X3*U_hat))) %*% matrix(c(X1*U_hat))
Omega_hat_PNW[1,2] = t(matrix(c(X1*U_hat))) %*% matrix(c(X2*U_hat))
Omega_hat_PNW[2,2] = t(matrix(c(X2*U_hat))) %*% matrix(c(X2*U_hat))
Omega_hat_PNW[3,2] = t(matrix(c(X3*U_hat))) %*% matrix(c(X2*U_hat))
Omega_hat_PNW[1,3] = t(matrix(c(X1*U_hat))) %*% matrix(c(X3*U_hat))
Omega_hat_PNW[2,3] = t(matrix(c(X2*U_hat))) %*% matrix(c(X3*U_hat))
Omega_hat_PNW[3,3] = t(matrix(c(X3*U_hat))) %*% matrix(c(X3*U_hat))
for( j in 1:m ){
 temp_Omega_hat_PNW = matrix(0,3,3)
 w = 1 - (j/(m+1))
 for( t in (j+1):T ){
  temp_Omega_hat_PNW[1,1] = temp_Omega_hat_PNW[1,1] + sum( diag( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[2,1] = temp_Omega_hat_PNW[2,1] + sum( diag( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[3,1] = temp_Omega_hat_PNW[3,1] + sum( diag( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[1,2] = temp_Omega_hat_PNW[1,2] + sum( diag( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[2,2] = temp_Omega_hat_PNW[2,2] + sum( diag( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[3,2] = temp_Omega_hat_PNW[3,2] + sum( diag( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[1,3] = temp_Omega_hat_PNW[1,3] + sum( diag( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[2,3] = temp_Omega_hat_PNW[2,3] + sum( diag( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) ) )
  temp_Omega_hat_PNW[3,3] = temp_Omega_hat_PNW[3,3] + sum( diag( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) ) )
 }
 Omega_hat_PNW = Omega_hat_PNW + 2 * w * temp_Omega_hat_PNW
}
Omega_hat_PNW = 1/(N*T) * Omega_hat_PNW
########################################################################
# ESTIMATE PANEL NEWEY-WEST OMEGA (ENDED)
########################################################################



########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################

### OMEGA HAT 1ST TERM #################################################
Omega_hat_1 = matrix(0,3,3)
for( i in 1:N ){
 Omega_hat_1[1,1] = Omega_hat_1[1,1] + sum( t(matrix(X1[i,]*U_hat[i,],1,T)) %*% matrix(X1[i,]*U_hat[i,],1,T) )
 Omega_hat_1[2,1] = Omega_hat_1[2,1] + sum( t(matrix(X2[i,]*U_hat[i,],1,T)) %*% matrix(X1[i,]*U_hat[i,],1,T) )
 Omega_hat_1[3,1] = Omega_hat_1[3,1] + sum( t(matrix(X3[i,]*U_hat[i,],1,T)) %*% matrix(X1[i,]*U_hat[i,],1,T) )
 Omega_hat_1[1,2] = Omega_hat_1[1,2] + sum( t(matrix(X1[i,]*U_hat[i,],1,T)) %*% matrix(X2[i,]*U_hat[i,],1,T) )
 Omega_hat_1[2,2] = Omega_hat_1[2,2] + sum( t(matrix(X2[i,]*U_hat[i,],1,T)) %*% matrix(X2[i,]*U_hat[i,],1,T) )
 Omega_hat_1[3,2] = Omega_hat_1[3,2] + sum( t(matrix(X3[i,]*U_hat[i,],1,T)) %*% matrix(X2[i,]*U_hat[i,],1,T) )
 Omega_hat_1[1,3] = Omega_hat_1[1,3] + sum( t(matrix(X1[i,]*U_hat[i,],1,T)) %*% matrix(X3[i,]*U_hat[i,],1,T) )
 Omega_hat_1[2,3] = Omega_hat_1[2,3] + sum( t(matrix(X2[i,]*U_hat[i,],1,T)) %*% matrix(X3[i,]*U_hat[i,],1,T) )
 Omega_hat_1[3,3] = Omega_hat_1[3,3] + sum( t(matrix(X3[i,]*U_hat[i,],1,T)) %*% matrix(X3[i,]*U_hat[i,],1,T) )
}
Omega_hat_1 = min(c(N,T))/N/(N*T^2) * Omega_hat_1

### OMEGA HAT 2ND TERM #################################################
Omega_hat_2 = matrix(0,3,3)
for( t in 1:T ){
 Omega_hat_2[1,1] = Omega_hat_2[1,1] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[2,1] = Omega_hat_2[2,1] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[3,1] = Omega_hat_2[3,1] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[1,2] = Omega_hat_2[1,2] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[2,2] = Omega_hat_2[2,2] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[3,2] = Omega_hat_2[3,2] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[1,3] = Omega_hat_2[1,3] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[2,3] = Omega_hat_2[2,3] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
 Omega_hat_2[3,3] = Omega_hat_2[3,3] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
}
Omega_hat_2 = min(c(N,T))/T/(N^2*T) * Omega_hat_2

### OMEGA HAT 3RD TERM #################################################
Omega_hat_3 = matrix(0,3,3)
Omega_hat_3[1,1] = t(matrix(c(X1*U_hat),,1)) %*% matrix(X1*c(1*U_hat),,1)
Omega_hat_3[2,1] = t(matrix(c(X2*U_hat),,1)) %*% matrix(X1*c(1*U_hat),,1)
Omega_hat_3[3,1] = t(matrix(c(X3*U_hat),,1)) %*% matrix(X1*c(1*U_hat),,1)
Omega_hat_3[1,2] = t(matrix(c(X1*U_hat),,1)) %*% matrix(X2*c(1*U_hat),,1)
Omega_hat_3[2,2] = t(matrix(c(X2*U_hat),,1)) %*% matrix(X2*c(1*U_hat),,1)
Omega_hat_3[3,2] = t(matrix(c(X3*U_hat),,1)) %*% matrix(X2*c(1*U_hat),,1)
Omega_hat_3[1,3] = t(matrix(c(X1*U_hat),,1)) %*% matrix(X3*c(1*U_hat),,1)
Omega_hat_3[2,3] = t(matrix(c(X2*U_hat),,1)) %*% matrix(X3*c(1*U_hat),,1)
Omega_hat_3[3,3] = t(matrix(c(X3*U_hat),,1)) %*% matrix(X3*c(1*U_hat),,1)
Omega_hat_3 = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3

### OMEGA HAT 4TH TERM FOR THOMPSON ####################################
m = 2
Omega_hat_4_thompson = matrix(0,3,3)
for( j in 1:m ){
 temp_Omega_hat_4_thompson = matrix(0,3,3)
 w = 1 - (j/(m+1))
 for( t in (j+1):T ){
  temp_Omega_hat_4_thompson[1,1] = temp_Omega_hat_4_thompson[1,1] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,1] = temp_Omega_hat_4_thompson[2,1] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[3,1] = temp_Omega_hat_4_thompson[3,1] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[1,2] = temp_Omega_hat_4_thompson[1,2] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,2] = temp_Omega_hat_4_thompson[2,2] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[3,2] = temp_Omega_hat_4_thompson[3,2] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[1,3] = temp_Omega_hat_4_thompson[1,3] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[2,3] = temp_Omega_hat_4_thompson[2,3] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_4_thompson[3,3] = temp_Omega_hat_4_thompson[3,3] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_4_thompson = 1 / (N^2*T) * temp_Omega_hat_4_thompson
 Omega_hat_4_thompson = Omega_hat_4_thompson + temp_Omega_hat_4_thompson
}
Omega_hat_4_thompson = min(c(N,T))/T * Omega_hat_4_thompson

### OMEGA HAT 5TH TERM FOR THOMPSON ####################################
m = 2
Omega_hat_5_thompson = matrix(0,3,3)
for( j in 1:m ){
 temp_Omega_hat_5_thompson = matrix(0,3,3)
 w = 1 - (j/(m+1))
 for( t in (j+1):T ){
  temp_Omega_hat_5_thompson[1,1] = temp_Omega_hat_5_thompson[1,1] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,1] = temp_Omega_hat_5_thompson[2,1] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[3,1] = temp_Omega_hat_5_thompson[3,1] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[1,2] = temp_Omega_hat_5_thompson[1,2] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,2] = temp_Omega_hat_5_thompson[2,2] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[3,2] = temp_Omega_hat_5_thompson[3,2] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[1,3] = temp_Omega_hat_5_thompson[1,3] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[2,3] = temp_Omega_hat_5_thompson[2,3] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_5_thompson[3,3] = temp_Omega_hat_5_thompson[3,3] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
 }
 temp_Omega_hat_5_thompson = w / (N^2*T) * temp_Omega_hat_5_thompson
 Omega_hat_5_thompson = Omega_hat_5_thompson + temp_Omega_hat_5_thompson
}
Omega_hat_5_thompson = min(c(N,T))/T * Omega_hat_5_thompson

### OMEGA HAT 6TH TERM FOR THOMPSON ####################################
m = 2
Omega_hat_6_thompson = matrix(0,3,3)
for( j in 1:m ){
 temp_Omega_hat_6_thompson = matrix(0,3,3)
 w = 1 - (j/(m+1))
 for( t in (j+1):T ){
  temp_Omega_hat_6_thompson[1,1] = temp_Omega_hat_6_thompson[1,1] + sum( matrix(X1[,t]*U_hat[,t],N,1) * (matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,1] = temp_Omega_hat_6_thompson[2,1] + sum( matrix(X2[,t]*U_hat[,t],N,1) * (matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[3,1] = temp_Omega_hat_6_thompson[3,1] + sum( matrix(X3[,t]*U_hat[,t],N,1) * (matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[1,2] = temp_Omega_hat_6_thompson[1,2] + sum( matrix(X1[,t]*U_hat[,t],N,1) * (matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,2] = temp_Omega_hat_6_thompson[2,2] + sum( matrix(X2[,t]*U_hat[,t],N,1) * (matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[3,2] = temp_Omega_hat_6_thompson[3,2] + sum( matrix(X3[,t]*U_hat[,t],N,1) * (matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[1,3] = temp_Omega_hat_6_thompson[1,3] + sum( matrix(X1[,t]*U_hat[,t],N,1) * (matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[2,3] = temp_Omega_hat_6_thompson[2,3] + sum( matrix(X2[,t]*U_hat[,t],N,1) * (matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
  temp_Omega_hat_6_thompson[3,3] = temp_Omega_hat_6_thompson[3,3] + sum( matrix(X3[,t]*U_hat[,t],N,1) * (matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
 }
 temp_Omega_hat_^_thompson = 1 / (N^2*T) * temp_Omega_hat_6_thompson
 Omega_hat_6_thompson = Omega_hat_6_thompson + temp_Omega_hat_6_thompson
}
Omega_hat_6_thompson = min(c(N,T))/T * Omega_hat_6_thompson

### OMEGA HAT 7TH TERM FOR THOMPSON ####################################
m = 2
Omega_hat_7_thompson = matrix(0,3,3)
for( j in 1:m ){
 temp_Omega_hat_7_thompson = matrix(0,3,3)
 w = 1 - (j/(m+1))
 for( t in (j+1):T ){
  temp_Omega_hat_7_thompson[1,1] = temp_Omega_hat_7_thompson[1,1] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) * (matrix(X1[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,1] = temp_Omega_hat_7_thompson[2,1] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) * (matrix(X1[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[3,1] = temp_Omega_hat_7_thompson[3,1] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) * (matrix(X1[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[1,2] = temp_Omega_hat_7_thompson[1,2] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) * (matrix(X2[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,2] = temp_Omega_hat_7_thompson[2,2] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) * (matrix(X2[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[3,2] = temp_Omega_hat_7_thompson[3,2] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) * (matrix(X2[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[1,3] = temp_Omega_hat_7_thompson[1,3] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) * (matrix(X3[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[2,3] = temp_Omega_hat_7_thompson[2,3] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) * (matrix(X3[,t]*U_hat[,t],N,1)) )
  temp_Omega_hat_7_thompson[3,3] = temp_Omega_hat_7_thompson[3,3] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) * (matrix(X3[,t]*U_hat[,t],N,1)) )
 }
 temp_Omega_hat_7_thompson = w / (N^2*T) * temp_Omega_hat_7_thompson
 Omega_hat_7_thompson = Omega_hat_7_thompson + temp_Omega_hat_7_thompson
}
Omega_hat_7_thompson = min(c(N,T))/T * Omega_hat_7_thompson

### COMPUTE M HAT ######################################################
S1 = colSums(X1*U_hat)[2:T]
S1lag = colSums(X1*U_hat)[1:(T-1)]
rho1_hat = cov(S1,S1lag)/var(S1lag)
S2 = colSums(X2*U_hat)[2:T]
S2lag = colSums(X2*U_hat)[1:(T-1)]
rho2_hat = cov(S2,S2lag)/var(S2lag)
S3 = colSums(X3*U_hat)[2:T]
S3lag = colSums(X3*U_hat)[1:(T-1)]
rho3_hat = cov(S3,S3lag)/var(S3lag)
m = (1.8171 * ((rho1_hat^2/(1-rho1_hat)^4+rho2_hat^2/(1-rho2_hat)^4+rho3_hat^2/(1-rho3_hat)^4)/((1-rho1_hat^2)^2/(1-rho1_hat)^4+(1-rho2_hat^2)^2/(1-rho2_hat)^4+(1-rho3_hat^2)^2/(1-rho3_hat)^4))^(1/3) * T^(1/3))
m = min(c(m,T-1))

### OMEGA HAT 4TH TERM #################################################
Omega_hat_4 = matrix(0,3,3)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_4 = matrix(0,3,3)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_4[1,1] = temp_Omega_hat_4[1,1] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,1] = temp_Omega_hat_4[2,1] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[3,1] = temp_Omega_hat_4[3,1] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X1[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[1,2] = temp_Omega_hat_4[1,2] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,2] = temp_Omega_hat_4[2,2] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[3,2] = temp_Omega_hat_4[3,2] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X2[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[1,3] = temp_Omega_hat_4[1,3] + sum( matrix(X1[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[2,3] = temp_Omega_hat_4[2,3] + sum( matrix(X2[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
   temp_Omega_hat_4[3,3] = temp_Omega_hat_4[3,3] + sum( matrix(X3[,t]*U_hat[,t],N,1) %*% t(matrix(X3[,t-j]*U_hat[,t-j],N,1)) )
  }
  temp_Omega_hat_4 = w / (N^2*T) * temp_Omega_hat_4
  Omega_hat_4 = Omega_hat_4 + temp_Omega_hat_4
 }
 Omega_hat_4 = min(c(N,T))/T * Omega_hat_4
}

### OMEGA HAT 5TH TERM #################################################
Omega_hat_5 = matrix(0,3,3)
if( m >= 1 ){
 for( j in 1:m ){
  temp_Omega_hat_5 = matrix(0,3,3)
  w = 1 - (j/(m+1))
  for( t in (j+1):T ){
   temp_Omega_hat_5[1,1] = temp_Omega_hat_5[1,1] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,1] = temp_Omega_hat_5[2,1] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[3,1] = temp_Omega_hat_5[3,1] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X1[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[1,2] = temp_Omega_hat_5[1,2] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,2] = temp_Omega_hat_5[2,2] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[3,2] = temp_Omega_hat_5[3,2] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X2[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[1,3] = temp_Omega_hat_5[1,3] + sum( matrix(X1[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[2,3] = temp_Omega_hat_5[2,3] + sum( matrix(X2[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
   temp_Omega_hat_5[3,3] = temp_Omega_hat_5[3,3] + sum( matrix(X3[,t-j]*U_hat[,t-j],N,1) %*% t(matrix(X3[,t]*U_hat[,t],N,1)) )
  }
  temp_Omega_hat_5 = w / (N^2*T) * temp_Omega_hat_5
  Omega_hat_5 = Omega_hat_5 + temp_Omega_hat_5
 }
 Omega_hat_5 = min(c(N,T))/T * Omega_hat_5
}

########################################################################
# ESTIMATE CHS OMEGA (BEGIN)
########################################################################



########################################################################
# VARIANCE ESTIMATION (BEGIN)
########################################################################

### DESIGN MATRIX Q HAT ################################################
Q_hat = t(cbind(matrix(c(X1),,1),matrix(c(X2),,1),matrix(c(X3),,1))) %*% cbind(matrix(c(X1),,1),matrix(c(X2),,1),matrix(c(X3),,1)) / (N*T)

### STANDARD ERRORS ####################################################
SE_IID = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_iid )                                                                                                                       )%*%solve(Q_hat)/(N*T) ) )
SE_PNW = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_PNW )                                                                                                                       )%*%solve(Q_hat)/(N*T) ) )
SE_CRi = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CRt = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_2 )                                                                                                                         )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CGM = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 )                                                                                             )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_THO = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4_thompson + Omega_hat_5_thompson - Omega_hat_6_thompson - Omega_hat_7_thompson ) )%*%solve(Q_hat)/min(c(N,T)) ) )
SE_CHS = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1 + Omega_hat_2 - Omega_hat_3 + Omega_hat_4 + Omega_hat_5 )                                                                 )%*%solve(Q_hat)/min(c(N,T)) ) )
########################################################################
# VARIANCE ESTIMATION (ENDED)
########################################################################



# MNW BOOTSTRAP ITERATIONS #############################################
num_boot_MNW = 1000
########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (BEGIN)
########################################################################
MNW_t_list = NULL
for( boot_iterations in 1:num_boot_MNW ){
rademacher = matrix(-1+2*(runif(N*T)<0.5),N,T)
U_hat_wild = rademacher * U_hat
Y_boot = beta_hat[1,1] * X1 + beta_hat[2,1] * X2  + beta_hat[3,1] * X3 + U_hat_wild
beta_hat_boot = solve( t(cbind(matrix(c(X1),,1),matrix(c(X2),,1),matrix(c(X3),,1))) %*% cbind(matrix(c(X1),,1),matrix(c(X2),,1),matrix(c(X3),,1)) ) %*% t(cbind(matrix(c(X1),,1),matrix(c(X2),,1),matrix(c(X3),,1))) %*% matrix(c(Y_boot),,1)
U_hat_boot = Y_boot - beta_hat_boot[1,1] * X1 - beta_hat_boot[2,1] * X2 - beta_hat_boot[3,1] * X3

### OMEGA HAT 1ST TERM #################################################
Omega_hat_1_boot = matrix(0,3,3)
for( i in 1:N ){
 Omega_hat_1_boot[1,1] = Omega_hat_1_boot[1,1] + sum( t(matrix(X1[i,]*U_hat_boot[i,],1,T)) %*% matrix(X1[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[2,1] = Omega_hat_1_boot[2,1] + sum( t(matrix(X2[i,]*U_hat_boot[i,],1,T)) %*% matrix(X1[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[3,1] = Omega_hat_1_boot[3,1] + sum( t(matrix(X3[i,]*U_hat_boot[i,],1,T)) %*% matrix(X1[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[1,2] = Omega_hat_1_boot[1,2] + sum( t(matrix(X1[i,]*U_hat_boot[i,],1,T)) %*% matrix(X2[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[2,2] = Omega_hat_1_boot[2,2] + sum( t(matrix(X2[i,]*U_hat_boot[i,],1,T)) %*% matrix(X2[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[3,2] = Omega_hat_1_boot[3,2] + sum( t(matrix(X3[i,]*U_hat_boot[i,],1,T)) %*% matrix(X2[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[1,3] = Omega_hat_1_boot[1,3] + sum( t(matrix(X1[i,]*U_hat_boot[i,],1,T)) %*% matrix(X3[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[2,3] = Omega_hat_1_boot[2,3] + sum( t(matrix(X2[i,]*U_hat_boot[i,],1,T)) %*% matrix(X3[i,]*U_hat_boot[i,],1,T) )
 Omega_hat_1_boot[3,3] = Omega_hat_1_boot[3,3] + sum( t(matrix(X3[i,]*U_hat_boot[i,],1,T)) %*% matrix(X3[i,]*U_hat_boot[i,],1,T) )
}
Omega_hat_1_boot = min(c(N,T))/N/(N*T^2) * Omega_hat_1_boot

### OMEGA HAT 2ND TERM #################################################
Omega_hat_2_boot = matrix(0,3,3)
for( t in 1:T ){
 Omega_hat_2_boot[1,1] = Omega_hat_2_boot[1,1] + sum( matrix(X1[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X1[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[2,1] = Omega_hat_2_boot[2,1] + sum( matrix(X2[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X1[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[3,1] = Omega_hat_2_boot[3,1] + sum( matrix(X3[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X1[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[1,2] = Omega_hat_2_boot[1,2] + sum( matrix(X1[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X2[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[2,2] = Omega_hat_2_boot[2,2] + sum( matrix(X2[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X2[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[3,2] = Omega_hat_2_boot[3,2] + sum( matrix(X3[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X2[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[1,3] = Omega_hat_2_boot[1,3] + sum( matrix(X1[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X3[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[2,3] = Omega_hat_2_boot[2,3] + sum( matrix(X2[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X3[,t]*U_hat_boot[,t],N,1)) )
 Omega_hat_2_boot[3,3] = Omega_hat_2_boot[3,3] + sum( matrix(X3[,t]*U_hat_boot[,t],N,1) %*% t(matrix(X3[,t]*U_hat_boot[,t],N,1)) )
}
Omega_hat_2_boot = min(c(N,T))/T/(N^2*T) * Omega_hat_2_boot

### OMEGA HAT 3RD TERM #################################################
Omega_hat_3_boot = matrix(0,3,3)

Omega_hat_3_boot[1,1] = t(matrix(c(X1*U_hat_boot),,1)) %*% matrix(X1*c(1*U_hat_boot),,1)
Omega_hat_3_boot[2,1] = t(matrix(c(X2*U_hat_boot),,1)) %*% matrix(X1*c(1*U_hat_boot),,1)
Omega_hat_3_boot[3,1] = t(matrix(c(X3*U_hat_boot),,1)) %*% matrix(X1*c(1*U_hat_boot),,1)
Omega_hat_3_boot[1,2] = t(matrix(c(X1*U_hat_boot),,1)) %*% matrix(X2*c(1*U_hat_boot),,1)
Omega_hat_3_boot[2,2] = t(matrix(c(X2*U_hat_boot),,1)) %*% matrix(X2*c(1*U_hat_boot),,1)
Omega_hat_3_boot[3,2] = t(matrix(c(X3*U_hat_boot),,1)) %*% matrix(X2*c(1*U_hat_boot),,1)
Omega_hat_3_boot[1,3] = t(matrix(c(X1*U_hat_boot),,1)) %*% matrix(X3*c(1*U_hat_boot),,1)
Omega_hat_3_boot[2,3] = t(matrix(c(X2*U_hat_boot),,1)) %*% matrix(X3*c(1*U_hat_boot),,1)
Omega_hat_3_boot[3,3] = t(matrix(c(X3*U_hat_boot),,1)) %*% matrix(X3*c(1*U_hat_boot),,1)
Omega_hat_3_boot = min(c(N,T))/(N*T)/(N*T) * Omega_hat_3

SE_CGM_boot = sqrt( diag( solve(Q_hat)%*%( EV_Correction( Omega_hat_1_boot + Omega_hat_2_boot - Omega_hat_3_boot ) )%*%solve(Q_hat)/min(c(N,T)) ) )
t_wild_boot = (beta_hat_boot - beta_hat ) / matrix(SE_CGM_boot,3,1)
MNW_t_list = cbind( MNW_t_list, t_wild_boot )
}
crit_MNW = apply(abs(MNW_t_list),1,quantile,p=0.95)
########################################################################
# WILD BOOTSTRAP CRITICAL VALUE (MNW2021) (ENDED)
########################################################################




RESULTS = cbind(beta_hat,SE_IID,SE_CRi,SE_CRt,SE_CGM,SE_MEN,SE_THO,SE_CHS)
rownames(RESULTS) = c("MKT","SMB","HML")
t_ratios = matrix(RESULTS[,1],3,7)/RESULTS[,2:8]


paste(calendar(first_period)," - ",calendar(last_period),"  N=",N," T=",T,sep="")
RESULTS
t_ratios